## ##

rm(list = ls())

##

library(tidyverse)
library(lfe)
library(stargazer)
library(lubridate)
library(broom)
source('code/function_tidy_felm.R')

## Def helper function to filter string vectors

str_filter <- function (string, pattern, neg = F) 
{
  if (neg == F) {
    string %>% .[str_detect(., pattern)]
  }
  else {
    string %>% .[!str_detect(., pattern)]
  }
}

## Get data 

civ <- read_rds("data/civey_party_replication.rds")

## Get weights 

weights_party <- read_rds('data/weights_civey_party.rds') %>% 
  mutate(month_10 = month_10 + 7)
weights_MIP <- read_rds('data/weights_civey_mip.rds') %>% 
  mutate(month_10 = month_10 + 7)

## Monthly intervals

start <- as.Date("2021-01-01")
end <- as.Date("2021-10-01")

dates_m <- floor_date(seq(start, end, by = 'month'), 
                      unit = "month") + 11

## Month indicators

civ <- civ %>% 
  mutate(month_10 = cut(date, dates_m)) %>% 
  mutate(month_10 = as.numeric(as.factor(month_10)) - 7) %>% 
  filter(!is.na(month_10))

## Create month dummies and interact with treatment 

civ_did <- civ %>%
  mutate(month_10 = month_10 + 7) %>% 
  mutate(month_dummy = month_10) %>% 
  fastDummies::dummy_cols(select_columns = 'month_dummy') %>% 
  dplyr::select(-`month_dummy_6`, -month_dummy) %>% 
  mutate_at(vars(matches("month_dummy")), list("weakly" = ~.*treat_weakly,
                                               "heavy" = ~.*treat_heavy)) %>% 
  dplyr::select(id, month_10, matches("weakly|heavy"), green, ags,
                educ, age, gender, state)  %>% 
  dplyr::select(-treat_weakly, -treat_heavy)

## Add weights 

civ_did <- civ_did %>% 
  left_join(weights_party)

## Indicator for panel respondents 
## Ie respondents that are observed in all periods

civ_did <- civ_did %>% 
  arrange(id, month_10) %>% 
  group_by(id) %>% 
  mutate(n_months = length(unique(month_10))) %>% 
  ungroup()

## Formula for regression 

f1 <- paste0("green ~ ",
             colnames(civ_did) %>% 
               str_filter("weakly") %>% paste0(collapse = '+'),
             '+',
             colnames(civ_did) %>% 
               str_filter("heavy") %>% paste0(collapse = '+'),
             "|ags + month_10 + age + educ + gender |0|id + ags") %>% 
  as.formula()

## Estimate

m1_all <- felm(f1, data = civ_did)
m1_only_panel <- felm(f1, data = civ_did %>% 
                        filter(n_months == 9))
m1_weights <- felm(f1, data = civ_did %>% 
                     filter(!is.na(w_mz)), 
                   weights = civ_did$w_mz[!is.na(civ_did$w_mz)])

## Extract quantities

m1_all <- m1_all %>% tidy_felm()
m1_only_panel <- m1_only_panel %>% tidy_felm()
m1_weights <- m1_weights %>% tidy_felm()

## Issue salience

civ_mip <- read_rds("data/civey_issue_replication.rds")

## Monthly intervals 

start <- as.Date("2021-01-01")
end <- as.Date("2021-10-01")
dates_m <- floor_date(seq(start, end, by = 'month'), 
                      unit = "month") + 11

## Month indicators 

civ_mip <- civ_mip %>% 
  mutate(month_10 = cut(date, dates_m)) %>% 
  mutate(month_10 = as.numeric(as.factor(month_10)) - 7) %>% 
  filter(!is.na(month_10))

## Create month dummies and interact with treatment 

civ_did <- civ_mip %>%
  mutate(month_10 = month_10 + 7) %>% 
  mutate(month_dummy = month_10) %>% 
  fastDummies::dummy_cols(select_columns = 'month_dummy') %>% 
  dplyr::select(-`month_dummy_6`, -month_dummy) %>% 
  mutate_at(vars(matches("month_dummy")), list("weakly" = ~.*treat_weakly,
                                               "heavy" = ~.*treat_heavy)) %>% 
  dplyr::select(id, month_10, matches("weakly|heavy"), climate, ags, 
                educ, age, gender, state) %>% 
  dplyr::select(-treat_weakly, -treat_heavy)

## Add weights 

civ_did <- civ_did %>% 
  left_join(weights_party)

## Indicator for panel respondents 
## Ie respondents that are observed in all periods

civ_did <- civ_did %>% 
  arrange(id, month_10) %>% 
  group_by(id) %>% 
  mutate(n_months = length(unique(month_10))) %>% 
  ungroup()

## Formula for regression 

f2 <- paste0("climate ~ ",
             colnames(civ_did) %>% 
               str_filter("weakly") %>% 
               paste0(collapse = '+'),
             '+',
             colnames(civ_did) %>% str_filter("heavy") 
             %>% 
               paste0(collapse = '+'),
             "|ags + month_10 + age + educ + gender|0|id + ags") %>% 
  as.formula()

## Estimate

m2_mip <- felm(f2, data = civ_did)
m2_mip_only_panel <- felm(f2, data = civ_did %>% 
                            filter(n_months == 9))
m2_mip_weights <- felm(f2, data = civ_did %>% 
                         filter(!is.na(w_mz)), 
                       weights = civ_did$w_mz[!is.na(civ_did$w_mz)])

## Get quantities

m2_mip <- m2_mip %>% tidy_felm()
m2_mip_only_panel <- m2_mip_only_panel %>% tidy_felm()
m2_mip_weights <- m2_mip_weights %>% tidy_felm()

## Plot

pd <- position_dodge(0.6)

#### Figure 4 - MIP and Green Party Support ####

m2_mip %>% 
  mutate(what = 'Climate change salience') %>% 
  bind_rows(m1_all %>% mutate(what = "Green party support")) %>% 
  filter(str_detect(term, "heavy|weakly")) %>% 
  mutate(period = substr(term, 13, 13),
         treat = substr(term, 15, 200)) %>% 
  mutate(period = as.numeric(period) - 7) %>% 
  ggplot(aes(period, estimate * 100, treat)) +
  geom_hline(yintercept = 0, linetype = 'dotted') +
  geom_vline(xintercept = -1, linetype = 'dotted') +
  geom_errorbar(aes(ymin = conf.low* 100, ymax = conf.high* 100, color = treat), 
                width = 0, position = pd)+
  geom_point(aes(color = treat, shape = treat), position = pd) +
  facet_wrap(~what) +
  theme_bw() +
  theme(legend.position = "bottom") +
  xlab("Months relative to flood") +
  ylab("Estimate (percentage points)") +
  scale_color_grey(end = 0.6, 
                   name = 'Flood intensity\n(vs. unaffected)',
                   labels = c('Highly affected', 'Weakly affected')) +
  scale_shape_manual(values = c(19, 15), 
                     name = 'Flood intensity\n(vs. unaffected)',
                     labels = c('Highly affected', 'Weakly affected')) +
  scale_x_continuous(breaks = -7:3)

#### Figure A.9 - respondents for which panel data is available ####

m2_mip_only_panel %>% 
  mutate(what = 'Climate change salience') %>% 
  bind_rows(m1_only_panel %>% mutate(what = "Green party support")) %>% 
  filter(str_detect(term, "heavy|weakly")) %>% 
  mutate(period = substr(term, 13, 13),
         treat = substr(term, 15, 200)) %>% 
  mutate(period = as.numeric(period) - 7) %>% 
  ggplot(aes(period, estimate * 100, treat)) +
  geom_hline(yintercept = 0, linetype = 'dotted') +
  geom_vline(xintercept = -1, linetype = 'dotted') +
  geom_errorbar(aes(ymin = conf.low* 100, ymax = conf.high* 100, color = treat), 
                width = 0, position = pd)+
  geom_point(aes(color = treat, shape = treat), position = pd) +
  facet_wrap(~what) +
  theme_bw() +
  theme(legend.position = "bottom") +
  xlab("Months relative to flood") +
  ylab("Estimate (percentage points)") +
  scale_color_grey(end = 0.6, 
                   name = 'Flood intensity\n(vs. unaffected)',
                   labels = c('Highly affected', 'Weakly affected')) +
  scale_shape_manual(values = c(19, 15), 
                     name = 'Flood intensity\n(vs. unaffected)',
                     labels = c('Highly affected', 'Weakly affected')) +
  scale_x_continuous(breaks = -7:3)

#### Figure A.11 - results using weights ####

m2_mip_weights %>% 
  mutate(what = 'Climate change salience') %>% 
  bind_rows(m1_weights %>% mutate(what = "Green party support")) %>% 
  filter(str_detect(term, "heavy|weakly")) %>% 
  mutate(period = substr(term, 13, 13),
         treat = substr(term, 15, 200)) %>% 
  mutate(period = as.numeric(period) - 7) %>% 
  ggplot(aes(period, estimate * 100, treat)) +
  geom_hline(yintercept = 0, linetype = 'dotted') +
  geom_vline(xintercept = -1, linetype = 'dotted') +
  geom_errorbar(aes(ymin = conf.low* 100, ymax = conf.high* 100, color = treat), 
                width = 0, position = pd)+
  geom_point(aes(color = treat, shape = treat), position = pd) +
  facet_wrap(~what) +
  theme_bw() +
  theme(legend.position = "bottom") +
  xlab("Months relative to flood") +
  ylab("Estimate (percentage points)") +
  scale_color_grey(end = 0.6, 
                   name = 'Flood intensity\n(vs. unaffected)',
                   labels = c('Highly affected', 'Weakly affected')) +
  scale_shape_manual(values = c(19, 15), 
                     name = 'Flood intensity\n(vs. unaffected)',
                     labels = c('Highly affected', 'Weakly affected')) +
  scale_x_continuous(breaks = -7:3)

